#include "stdafx.h"
#include "PWM.h"

//Static const string
const char* CPWM::NT = "ACGT";
CPWM::CPWM(void)
{
    this->initialiation(0);
}
CPWM::CPWM(unsigned int iKmerLength)
{
    this->initialiation(iKmerLength);

}
CPWM::~CPWM(void)
{
    for (unsigned int i = 0; i < this->NUM_NT; i++) {
        delete [] this->caQualityScore[i];
    }
}
int CPWM::initialiation(unsigned int iKmerLength)
{
    if (iKmerLength > 0) {
        this->iKmerLength = iKmerLength;
        for (unsigned int i = 0; i < NUM_NT; i++) {
            caQualityScore[i] = new char[iKmerLength];
        }
    } else {
        this->iKmerLength = 0;
        for (unsigned int i = 0; i < this->NUM_NT; i++) {
            this->caQualityScore[i] = NULL;
        }
    }
    return(0);
}
// The version that read directly from file.
bool CPWM::readinQaulityScore(ifstream &ifile)
{
    char buffer[MAX_LINE]; //read in a line of score for each kmer
    if (ifile.good()) {
        ifile.getline(buffer, MAX_LINE);
        return(getQaulityScores(buffer));
    }
    return(false); //fail to get score
}

// Override the previous function using inputBuffer
bool CPWM::readinQaulityScore(FileInputBuffer* inputBuffer)
{
    char buffer[MAX_LINE]; //read in a line of score for each kmer
    if (inputBuffer != NULL) {
        inputBuffer->Getline(buffer, MAX_LINE);
        return(getQaulityScores(buffer));
    }
    return(false); //fail to get score
}

/**
* Get the most likely sequence from the quality score.
* The buffer it pass in is to store the returned sequence.
* Return the sum of Q score for the most likely happend sequence
*/
int CPWM::getSequence(char* caSeq)
{
    if (caSeq != NULL && this->caQualityScore != NULL) {
        unsigned int i;
        int scoreSum = 0;
        for (i = 0; i < this->iKmerLength; i++) {
            char maxScore = this->SCORE_LOWERBOUND;
            for (unsigned int j = 0; j < this->NUM_NT; j++) {
                if (maxScore < this->caQualityScore[j][i]) {
                    maxScore = this->caQualityScore[j][i];
                    caSeq[i] = this->NT[j];
                }
            }
            scoreSum += maxScore;
        }
        caSeq[i] = '\0';
        return(scoreSum);
    } else {
        return(this->SCORE_LOWERBOUND*this->iKmerLength); //A indicater that one may not use
    }
}
/**
* Get sum of scores of the most likely sequence of the read.
* This function is very similar to getSequence except record the most likely sequence.
*/
int CPWM::getBestScoreSum()
{

    int scoreSum = 0;
    for (unsigned int i = 0; i < this->iKmerLength; i++) {
        char maxScore = this->SCORE_LOWERBOUND;
        for (unsigned int j = 0; j < this->NUM_NT; j++) {
            if (maxScore < this->caQualityScore[j][i]) {
                maxScore = this->caQualityScore[j][i];
            }
        }
        scoreSum += maxScore;
    }
    return(scoreSum);
}


/** private function called by readinQaulityScore(), with
 *  ptr buffer should contain the score info in K teraplet form */
bool CPWM::getQaulityScores(char* buffer)
{
    char* pch;
    if (buffer != NULL) {
        int score = 0; //tmp int to store the score for each bases
        for (unsigned int i = 0; i < this->iKmerLength; i++) {
            for (unsigned int j = 0; j < this->NUM_NT; j++) {
                if (i == 0 && j == 0) {
                    pch = strtok(buffer, " \t");
                } else {
                    pch = strtok(NULL, " \t");
                }
                if (pch != NULL) {
                    score = atoi(pch); //Change the string to integer score
                } else {
                    cout << "End of reads file" << endl;
                    return(false); //fail to get score
                }
                if (score >= this->SCORE_LOWERBOUND && score <= this->SCORE_UPPERBOUND) {
                    this->caQualityScore[j][i] = (char)(score); //convert the integer score to a char
                } else {
                    cout << "Invalid score" << endl;
                    return(false); //fail to get score
                }
            }
        }
        return(true); // successfully get score
    }
    return(false); //fail to get score
}

/**
* function will calculate the posible scores and their frequency based on a PWM
* It reaturned a map with score and frequency. It should be sorted by score.
* Should be completed by adding up transition matrix
*/
map<int, double>* CPWM::getScoresHist(void)
{
    map<int, double>* scoreFreqMap = new map<int, double>;
    (*scoreFreqMap)[0] = 1.0; //initialization

    for (unsigned int i = 0; i <  this->iKmerLength; i++) {
        // getSumsOfScoresFrequencyInDPStep will delete the old map and return a new address
        scoreFreqMap = getSumsOfScoresFrequencyInDPStep(scoreFreqMap, i);
    }
    return scoreFreqMap;
}

// To do one step in Tade's DP algorithm in finding the combination of sum of scores
map<int, double>* CPWM::getSumsOfScoresFrequencyInDPStep(map<int, double>* ScoreMap, int steps)
{

    map<int, double>::reverse_iterator it;
    map<int, double>* NewScoreMap = new map<int, double>;
    unsigned int rankSize = 0;

    /* instead of using iterator from the begining (increasing), iterator go back ward (decreasing) */
    for (it = ScoreMap->rbegin(); it != ScoreMap->rend(); it++) {
        for (unsigned int i = 0; i < this->NUM_NT ; i++) {
            int record = it->first + (int)this->caQualityScore[i][steps]; //a instance of posible score
            if (NewScoreMap->find(record) == NewScoreMap->end()) { //check if it is already in hist
                if (record > 0 || NewScoreMap->size() < SIZE_OF_RANKED_SCORE_MAP) { /* Set a threshold for */
                    (*NewScoreMap)[record] = it->second;
                }
            } else {
                (*NewScoreMap)[record] += it->second;
            }
        }
        rankSize++;
        if (rankSize > SIZE_OF_RANKED_SCORE_MAP) //Limiate the rank in our constant
            break;
    }
    /* check
    for(it = ScoreMap->rbegin(); it != ScoreMap->rend(); it++) {
    	cout<<it->first<<","<<it->second<<endl;
    }cout<<endl; */
    delete ScoreMap;
    return NewScoreMap;
}

double CPWM::getSignificance(map<int, double>* hist, int score)
{

    map<int, double>::reverse_iterator it;

    double count = 0.0;
    /* instead of using iterator from the beginning (increasing),
     *  iterator go back ward (decreasing) */
    for (it = hist->rbegin(); it != hist->rend(); it++) {
        // int oneRecordScore = it->first;
        // double oneRecordFreq = it->second; //For debug
        if (score <= it->first) {
            count += (double)it->second;
        } else {
            break;
        }
    }
    return(count); //How many path has better score
}

double CPWM::getPValue(map<int, double>* hist, int score)
{
    double count = this->getSignificance(hist, score);
    double totalCount = pow((double)this->NUM_NT, (double)this->iKmerLength);
    return(count / totalCount);
}

/**
* A function use for test set up: Given a kmer set the correspondent NT SCORE_UPPERBOUND and other SCORE_LOWERBOUND
*/
void CPWM::setPWM(char* kmer)
{
    for (unsigned int i = 0; i < this->NUM_NT; i++) {
        for (unsigned int j = 0; j < this->iKmerLength; j++) {
            this->caQualityScore[i][j] = this->SCORE_LOWERBOUND;
        }
    }
    for (unsigned int i = 0; i < this->iKmerLength; i++) {
        this->caQualityScore[this->NTMap(kmer[i])][i] = this->SCORE_UPPERBOUND;
    }
}

//Set pwm According to position for test purpose
void  CPWM::setPWM(int iPosi, double qA, double qC, double qG, double qT)
{
    this->caQualityScore[0][iPosi] = convertProb2QScore(qA);
    this->caQualityScore[1][iPosi] = convertProb2QScore(qC);
    this->caQualityScore[2][iPosi] = convertProb2QScore(qG);
    this->caQualityScore[3][iPosi] = convertProb2QScore(qT);
}

unsigned int CPWM::NTMap(char base)
{
    switch (base) {
    case 'a':
    case 'A':
        return(0x00);
    case 'c':
    case 'C':
        return(0x01);
    case 'g':
    case 'G':
        return(0x02);
    case 't':
    case 'T':
        return(0x03);
    default:
        return(0x04);
    }
}
/* Return the quality score of a base and position
 * with 0, 1, 2, 3 as query */
char CPWM::getScore(unsigned int ntId, unsigned int posId)
{
    if (ntId < CPWM::NUM_NT)
        return (this->caQualityScore[ntId][posId]);
    else
        return (CPWM::SCORE_LOWERBOUND);
}
/* Return the quality score of a base and position
 * with 'A', 'C', 'G' or 'T' as query */
char  CPWM::getQScore(char nt, unsigned int posId)
{
    unsigned int ntId = this->NTMap(nt);
    if (ntId < CPWM::NUM_NT)
        return(caQualityScore[ntId][posId]);
    else
        return (CPWM::SCORE_LOWERBOUND);
}
/**
* Given the mapped refrence, return the sum of the quality score for the correpondent nucleotide
*/
int CPWM::getScoreSum(char* mappedRef)
{
    int ScoreSum = 0;
    for (unsigned int i = 0; i < this->iKmerLength && i < (unsigned int)strlen(mappedRef); i++) {
        unsigned int ntId = this->NTMap(mappedRef[i]);
        if (ntId < this->NUM_NT)
            ScoreSum += (int)this->caQualityScore[ntId][i];
        else
            ScoreSum += this->SCORE_LOWERBOUND;
    }
    return(ScoreSum);
}
double CPWM::convertQScore2Prob(char score)
{
    switch (score) {
    case 40:
        return 1.0;
    case -40:
        return 0;
    case -5:
        return 0.25;
    default:
        double a = (double)score / -10.0;
        double r = pow(10.0, a);
        double errorProb = r / (1 + r);
        return(1 - errorProb);
    }
}
char CPWM::convertProb2QScore(double Prob)
{
    double errorProb = 1 - Prob;
    double dScore = -10 * (log(errorProb / (1 - errorProb)) / log(10.0));
    int iScore = (int)floor(dScore);
    char cScore = (char)(iScore % 128);
    return(cScore);
}

int CPWM::getQualityScore(CPWM &pwm, char* caRef)
{
    int qualityScore = pwm.getScoreSum(caRef);
    if (qualityScore < 0) {
        ERR //Check if this is really a bad Read
        cout << " Sum of quality score less than one" << endl;
        for (unsigned int i = 0; i < pwm.iKmerLength; i++) {
            cout << '(';
            for (unsigned int j = 0; j < pwm.NUM_NT; j++) {
                cout << (int)pwm.getScore(j, i) << ',';
            }
            cout << ')' << endl;
        }
        cout << '\n' << endl;
    }
    return (qualityScore);
}


